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This paper considers the aeroelastic optimization of a subsonic transport wingbox under 
a variety of static and dynamic aeroelastic constraints. Three types of design variables are 
utilized: structural variables (skin thickness, stiffener details), the quasi-steady deflection 
scheduling of a series of control surfaces distributed along the trailing edge for maneuver 
load alleviation and trim attainment, and the design details of an LQR controller, which 
commands oscillatory hinge moments into those same control surfaces. Optimization 
problems are solved where a closed loop flutter constraint is forced to satisfy the required 
flight margin, and mass reduction benefits are realized by relaxing the open loop flutter 
requirements. 


I. Introduction 

A eroelastic tailoring of wing structures has typically consisted of optimizing the thickness distribution of the 
various wing box components (skins, spars, etc.), or optimizing the stacking sequence, in the case of composite 
structures [1]. Design metrics may include static aeroelastic trim stresses, panel buckling, aeroelastic stability 
(divergence, flutter), gust response, maneuver/trim drag, etc. In addition to structural sizing variables, one may also 
consider distributed multiple control surfaces (DMCS) along the trailing edge of the wing (and, potentially, the 
leading edge). Optimizing the static and dynamic deflections of these control surfaces would change the apparent 
stiffness and mass of the wing structure, rather than the true quantities, but these divisions are blurred in any case, 
due to the coupled nature of an aeroelastic system. Distributed control surfaces may be equally adept at altering 
load paths and dynamic behavior of a wing structure for improved performance in terms of any of the design metrics 
mentioned above. Furthermore, simultaneously designing the wing structure and the actuator deflections will allow 
an optimizer to take advantage of the synergies between the two types of design variables. 

Distributed multiple control surfaces may be used in two distinct ways [2]. First, their deflection patterns may 
be scheduled in a quasi-steady manner for maneuver load alleviation. A separate deflection pattern may be 
optimized for each trim load case, where the distributed control surfaces would simultaneously alleviate the loads 
into the wing structure and help maintain aerodynamic trim. This may be done by solving an over- determined trim 
problem, where flap deflection scheduling is optimized such that trim is maintained (posed as a constraint) and some 
objective of interest is minimized: root bending moment [3], deformations at given finite element nodes [4], drag, 
etc. This typically is posed as a nested optimization solution, where structural design variables are applied to the 
over-determined trim solution. Alternatively, structural design variables may be optimized concurrently with the 
control surface scheduling [5], either iteratively or simultaneously. This latter scenario is of interest for this work. 

Secondly, the control surfaces may be used for active flutter mitigation or gust rejection, where each flap 
dynamically oscillates about the steady-state position found with the methods in the previous paragraph 
(superposition). This introduces a 3 rd type of design variable into the aeroelastic optimization (in addition to the 
structural variables and steady flap deflection variables): the details of the controller. Though aeroservoelastic 
optimization may be accomplished through a variety of means [2], one convenient technique is to use a linear 
quadratic controller, and optimize the details of the weighting matrices in the quadratic performance index 
(controller cost function). This is demonstrated in Refs. [6], [7], and [8]. Furthermore, the performance index itself 
may be used as an objective function or constraint [9][10][1 1] during the optimization process. Papers that study the 
use of multiple and/or distributed control surfaces for aeroservoelastic control (though not for optimization) can be 
found in Refs. [12], [13], and [14]. 
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Aeroservoelasticity in general is concerned with flutter, gust alleviation, ride quality, etc. Of these dynamic 
performance metrics, the current paper only considers flutter, in an effort to understand the structural weight savings 
that may be gained when an open loop flutter margin constraint is relaxed, and closed loop control is instead used to 
impose the required safety margin. This work is conducted on the un-deflected Common Research Model (uCRM) 
[15], a jig shape of the flying wing model developed in Ref. [16]. The uCRM is a generic subsonic transport 
configuration with a standard rib-spar wingbox, and has been outfitted for this work with a series of control surfaces 
distributed along the trailing edge. 

As noted above, three types of design variables are considered here, which in this paper are introduced into the 
aeroelastic optimization problem in stages. First (Section IV), only structural design variables (skin thickness and 
stiffener geometry for the skins, ribs, and spars) are used to minimize the mass of the wingbox, under static 
aeroelastic constraints only. Next (Section V), the quasi-steady deflections of each control surface (DMCS), at each 
maneuver load case, are added to the vector of design variables, and the optimization problem is re-solved. The 
impact (on the minimum attainable mass) of including an open loop flutter constraint is then assessed (Section VI). 
Finally (Section VII), controller design variables are added to the design variable vector, and a closed loop flutter 
constraint is imposed, along with the open loop constraint. Trade-offs (Pareto fronts) between optimal wing mass, 
open loop flutter margin, and control cost (performance index) are all presented. 

II. The Common Research Model 

All of the work in this paper is conducted on the conceptual Common Research Model (CRM). The model 
presented in Ref. [16] is a lg flying shape suitable for rigid aerodynamic analysis. An un-deflected jig shape version 
(uCRM) developed in Ref. [15] is suitable for aeroelastic analysis, and is used here. This transonic transport 
configuration has a wing span of 58.7 m, a mean aerodynamic chord of 7.0 m, an aspect ratio of 9, a taper ratio of 
0.275, a sweep angle of 35°, and a cruise Mach number of 0.85. 

The topology of the wingbox developed in Ref. [15] is also used here, seen in Figure 1. This structure consists 
of an upper skin, a lower skin, a leading edge spar (located at 10% chord at the root and 35% at the tip), a trailing 
edge spar (60% at both the root and the tip), and 43 ribs oriented perpendicular to the leading edge. All shell 
members (ribs, spars, skins) are outfitted with blade stiffeners, where the pitch is equal to 0.33 m throughout. Run- 
out stiffeners are utilized down the span: 18 stiffeners are found for panels and ribs at the root, but only 3 at the tip. 
The stiffeners are not modeled explicitly, but instead smeared into the shell stiffness properties. The entire wing 
structure is constructed of aluminum ( E = 70 GPa, v = 0.3, p = 2780 kg/m 3 , <j Y = 330 MPa), and discretized into 
16,000 triangular finite elements. 



Figure 1. uCRM wingbox topology and outer mold line, taken from Ref. [15]. 

Twenty simple flap-like control surfaces are distributed along the trailing edge of the wing from root to chord, 
as shown in Figure 2. A hinge spring attached to each control surface is connected back to the trailing spar via 
interpolation elements. The mass of each DMCS is assumed proportional to the volume enclosed within the surface 
wedge, and the stiffness of each hinge spring is sized such that all of the control surfaces have the same natural 
vibration frequency. The structure within the leading edge of the wing is not explicitly modeled, though an inertial 
effect is captured with a series of lumped masses attached to the wing box via interpolation elements. These masses, 
along with similar representations for the engine and the fuel, are also shown in Figure 2. In addition to the inertial 
loads from the lumped masses, a thrust force is applied to the central engine node. 

Aerodynamic paneling for the wing, horizontal tail, vertical tail, fuselage, and engine (the latter two represented 
as cruciforms) is shown in Figure 3, with a total of 4,600 panels. For static aeroelastic trim analysis, the entire 
vehicle representation of Figure 3 is utilized, and trailing edge control surface hinge springs are modeled as rigid. In 
this case, commanded quasi-steady flap deflections (relative to the wing) are exactly obtained. For dynamic flutter 
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analysis, only the wing panels are utilized. Furthermore, the control surfaces are now allowed to vibrate freely, and 
commanded flap deflections may be resisted by aerodynamic and inertial hinge moments. 



Figure 2. Control surfaces attached to the uCRM model, and lumped mass representation of engine and fuel 

loading. 



Figure 3. Aerodynamic paneling of the uCRM. 
III. Static Aeroelastic Modeling 


A. Trim Loads 

The shell finite elements used to model the wing structure are defined by a combination of linear strain triangles 
(LST) and discrete Kirchhoff triangles (DKT) [17]. For static airloads, a linear vortex lattice method [18] is used to 
model the aerodynamic lifting surfaces. A finite plate spline (FPS) method [19] is used to transfer downwash and 
pressures between the aerodynamic and structural/control surface modules. Only information pertaining to the wing 
is transferred in this way: the remaining aerodynamic surfaces are not explicitly tied to any structure. 

The wing box structure is sized across two different types of static maneuvers. The first type is a longitudinal 
maneuver (pull-up, push-over), where the system is trimmed via the angle of attack a and the elevator deflection 6. 
The trimmed values of these two quantities are automatically found by augmenting the aero-structural coupling 
equations. The final equation for a longitudinal maneuver is written as: 
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The first row of Eq. 1 is the finite element analysis: K is the stiffness matrix, and the solution vector x s has six 
degrees of freedom per node (three displacements and three rotations). The s subscript indicates a “symmetric” 
term, to differentiate from the anti- symmetric terms below. Forcing functions include self-weight inertial loading 
Fg r av (scaled by the maneuver load factor N, and accounting for both the weight of the wing structure and the 
lumped masses in Figure 2), thrust loading F thrust from the engine, and aerodynamic forces. Aerodynamic forces 
are written as q ■ Q ■ C p , where C p is a vector of differential pressure coefficients acting on each panel of the vehicle, 
Q is an interpolation function derived from FPS, and q is the dynamic pressure. 

The second row of Eq. 1 is the aerodynamic analysis, where D s is the matrix of aerodynamic influence 
coefficients (AIC) and the subscript indicates a symmetric aerodynamic condition about the centerline of the 
airplane in Figure 3. This equation is driven by several terms: downwash due to angle of attack L a ■ a (where L a is 
a linear operator that converts the scalar angle of attack into a downwash at each panel), elevator deflection L s ■ 6, 
built in camber/twist of the wing and tail jig shapes z jig , and downwash induced by structural wing deformation. 
This latter term is written as P ■ x s , where P is a second interpolation function, also derived from FPS-based 
methods. A final downwash term is needed if distributed control surfaces are used along the trailing edge during the 
maneuver. The deflection of each control surface is grouped into the vector y s , and L y is a matrix that converts 
these deflections into the appropriate downwash at each panel. 

These deflections y s are known, specified quantities during the solution of Eq. 1 : a and 6 are found that trim 
the system in the presence of these and other terms. Trim equations are written in the 3 rd and 4 th rows of Eq. 1 : 
q ■ S L and q ■ S m convert the differential pressure vector C v into a total aerodynamic lift and aerodynamic pitching 
moment (about the aircraft center of gravity). Lift must offset the total weight of the vehicle ( N ■ W), and the net 
pitching moment must be zero. 

A second type of static maneuver considered here is a rolling trim analysis (Eq. 2), where the deflection /? of an 
outboard wing aileron is found such that a constant specified non-dimensional roll rate p ■ L/U is maintained, with 
no rolling acceleration. In this analysis, p is the dimensional roll rate, L is the semi-span, and U is the flight speed. 
The “aileron” is formed by grouping the 4 control surfaces placed between 70% and 90% of the semi-span in Figure 
2. The system is simultaneously trimmed longitudinally for steady level flight (N= 1) with the angle of attack a. 
The rolling analysis requires an anti-symmetric condition about the centerline of the airplane; the level flight 
anlaysis uses a symmetric condition. Owing to the linear nature of the methods used here, structural deformations 
and aerodynamic pressures for the two conditions can be solved separately, and then added together to obtain the 
total aeroelastic state. The final equation for this maneuver (rolling at level flight) is written as: 
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The a subscripts indicates an “anti-symmetric” term, q ■ S p converts the aerodynamic pressures into a rolling 
moment about the centerline (which is ultimately set to zero for a constant roll rate), and L p converts the rolling 
motion into a downwash term. It is also seen that distributed control surfaces along the trailing edge, if used, are 
separated into symmetric and anti- symmetric scheduling. Trailing edge control surfaces used to define the (3 
rotation in the trim module are excluded from the y a terms. 

B. Stresses 

Having solved the static trim Eqs. 1 and 2, stresses and strains can be computed for each static load case, and a 
knock-up (safety) factor of 30% is applied to each elemental stress value. The von Mises failure function is 
computed for each finite element. The spar, rib, and skin structures are then divided into patches (seen in Figure 1, 
with a total of 217 patches), and the Kreisselmeier-Steinhauser (KS) function [20] is used to compress all of the 
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elemental failure function values within a given patch into a single metric. If all of the stress values within a patch 
are within their failure envelope, the KS function for that patch will be less than one. 

C. Panel Buckling 

After the stress analysis, buckling analyses are run for each stiffened panel in the upper and lower skins. 
Following Refs. [21] and [22], each buckling analysis is conducted with a Rayleigh-Ritz method (assumed buckling 
modes). Both global buckling of a stiffened panel (bordered by ribs and spars) and local buckling in between each 
stiffener is computed, where simply supported boundary conditions are used for both scenarios. The buckling 
equation is given as: 

(K s + l i n -R)-v n = 0 (3) 

where v n is the eigenvector (vector of modal amplitudes) associated with the n th eigenvalue fi n . The panel stiffness 
matrix K is computed based upon the stiffness properties defined at the wing level. Smeared stiffness properties are 
necessarily used for global panel buckling, whereas un-smeared shell properties are used for local inter- stringer 
buckling. The panel’s geometric stiffness matrix K s is assembled with element stresses computed from Eqs. 1 or 2. 

An eigenvalue fj. n greater than one indicates that, for the trimmed aeroelastic state, the panel has buckled. Each 
buckling computation for a given stiffened panel is compressed into a single value using a Kreisselmeier- 
Steinhauser function. For a panel with 3 stringers (such as found at the wing tip), 4 local buckling computations are 
performed, along with 1 global buckling computation. Since three modes (n = 3) are computed for each 
eigenproblem, a total of 15 eigenvalues are compressed into a single KS function for this example. The total 
number of stiffened panels in Figure 1 along the upper and lower skins is 88, which is also equal to the number of 
KS buckling values. 

D. Design Variables and Analytical Sensitivities 

Structural design variables for static aeroelastic optimization include shell thickness, stiffener thickness, and 
stiffener height, where simple blade stiffeners are used. These variables will have an impact on the wing and panel 
stiffness matrices K , K , and K s , the inertial loading F grav , the vehicle weight VF, and the term S m (which is 
dependent upon the center of gravity). In the event that the scheduling of distributed control surface deflections are 
used during the optimization, the relevant design variables are y s and y a . One set of symmetric deflections y s is 
designed for each load case (Eq. 1 or 2), and a set of anti-symmetric deflections ( y a ) are further included for each 
rolling load case (Eq. 2). 

Derivatives of the static aeroelastic response (stress and buckling aggregation parameters for each load case) 
with respect to all design variables are computed with the adjoint method. 

IV. Static Aeroelastic Optimization with Structural Design Variables 

The optimal distribution of patch-based skin thickness, stiffener thickness, and stiffener height variables (using 
the rib, spar, and skin design patches in Figure 1) are found that minimizes the wing mass (based on the volume of 
the finite element model) subject to aeroelastic strength and buckling constraints: 

min 

mass 

R 

r 0<q<l ( 4 ) 

KS c i < 1 i = 1 (N a ■ N l ) 

l KS^<1 i = 1 (/V M • N l ) 

where q are the design variables that have been appropriately normalized between 0 and 1 . Dimensionally, shell 
thicknesses are allowed to range between 3 mm and 30 mm, stiffener thicknesses between 3 mm and 30 mm, and 
stiffener heights between 30 mm and 150 mm. 

N l is the number of static load cases, KS a are the stress aggregation parameters (N a per load case), KS g are the 
buckling aggregation parameters (N g per load case). Three static load cases are considered for all optimization 
problems in this paper. Each has a Mach number of 0.85 and a dynamic pressure of 13.8 kPa. The first two load 
cases are longitudinal maneuvers (Eq. 1), with load factors of 2.5 and -1, respectively. The third load case is 
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dictated by Eq. 2, where the desired non-dimensional roll rate is 0.065. All cases are run with half-fuel, and the 
distributed control surfaces along the trailing edge (except the aileron deflection (3 used for lateral trim) are ignored: 
Ys = Ya = ® Quasi-steady control surface deflection design variables will be considered in section V. 

The design variables q are passed through a linearly-decaying cone-shape filter [23] prior to conversion into 
structural quantities, in order to prevent the difference in stiffness between adjacent patches from being too large. 
All optimization problems are solved with the Globally-Convergent Method of Moving Asymptotes [24]. The 
optimal result is shown in Figure 4, in terms of the shell thickness and smeared stiffener height of each design patch. 
The smeared stiffener thickness of each patch resides at the prescribed lower bound of 3 mm, and is not shown. The 
optimizer places the majority of the material in the upper and lower skins, with thickness and stiffener height tapered 
towards the lower bound at the wing tip, due to the drop in bending moment. Much larger stiffeners are utilized in 
the upper skins as compared to the lower skins, which is a buckling-driven result (discussed below). Rib design 
variables are all pushed to the lower bound, as are many of the spar variables, with the exception of the thick trailing 
edge root spar. Crushing loads are not included here, which may drive the optimizer to strengthen the rib members. 



Figure 4. Optimal structural design variables for the static aeroelastic optimization problem. 

Structural deflections for each load case are given in Figure 5, element-level von Mises failure functions in 
Figure 6, and KS buckling parameters in Figure 7. It is noted that the data in Figure 6 are not explicitly provided to 
the optimizer, but first compressed into KS 0 functions for each design segment. The buckling data in Figure 7 is 
KS and therefore is provided to the optimizer in the form shown. As discussed above, each KS^ is an aggregation 
of each global and local (inter-stringer) buckling eigenvalue for a given stiffened panel. For both stress data and 
buckling data, a value greater than one would be an infeasible aeroelastic constraint, and is not seen in either plot. 

load case 1 : 2.5g 
load case 3: roll 


jig shape 
load case 2: - ig 


Figure 5. Structural deflections of the optimal design in Figure 4 across each load case. 

Wing deflections in Figure 5 follow the expected pattern, though the rolling load case 3 has deflections nearly 
as large as the pull-up maneuver (2.5g). Furthermore, this trimmed rolling maneuver does not qualitatively show the 
torsion-dominated response that may be expected. The stress distribution (Figure 6) is markedly different however, 
with peak skin stresses shifted further outboard than seen for the longitudinal maneuvers. Farge portions of the 
outboard spars are stressed for this rolling maneuver as well, whereas for the longitudinal maneuvers, high spar 
stresses are confined to the engine attachment points. For load case 2 (-lg push over), these engine stresses are the 
critical values for the entire load case: stresses in the skins are relatively low. For load case 1, however, the majority 
of the upper and lower skin stresses are close to the allowable limit, which is a nearly fully-stressed design. Rib 
stresses are uniformly low, again likely due to the exclusion of crushing loads. 
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stress failure index 


0 0.2 0.4 0.6 0.8 1 

Figure 6. Stress-based failure indices of the optimal design in Figure 4 across each load case. 

Panels with negative KS buckling eigenvalue functions are left blank in Figure 7, as this is a physically 
meaningless result: for a given maneuver case, the loads would have to be entirely reversed for these panels to 
buckle. This clearly happens through most of the lower surface during the positive load factors (cases 1 and 3), and 
through most of the upper surface during the negative load factor (case 2). Otherwise, most of the panels have a 
buckling KS factor nearly equal to unity (active constraints) for the upper and lower skins of load case 1 and 2, 
respectively, as well as the outboard portion of load case 3 (where the stresses are large). As such, the buckling 
constraints in Eq. 4 are equally- strong design drivers as compared to the stress constraints, if not stronger. These 
buckling constraints are also presumably the cause of the large stiffeners utilized in the upper skins (Figure 4), as the 
2.5g maneuver causes stronger compressive forces in the upper skin than the -lg maneuver does in the lower skin. 

The mass of the optimal structure in Figure 4 is 6,175 kg; this structure is based on the volume of the finite 
element model (including smeared stiffeners), and does not include the lumped mass or control surfaces seen in 
Figure 1 and Figure 2. This value can be compared to values obtained below, when additional design variables and 
aeroelastic constraints are added to Eq. 4. 


lower skin 


upper skin 


load case 1 : 2. 


load case 2: - 


load case 3: rol l 



Figure 7. Buckling failure indices (KS functions) for each stiffened panel of the optimal design in Figure 4 

across each load case. 


V. Quasi-Steady DMCS Deflection Design Variables 

This section augments the optimization problem of Eq. 4 to include trailing edge flap deflection design 
variables (y s and y a ), in addition to the structural design variables. Simultaneous handling of wing structure design 
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variables and control surface scheduling design variables, within the same optimization loop, is in contrast to most 
of the existing literature on this topic [5]. Control surface variables are typically computed within an inner loop 
(over-determined trim optimization), and structural design variables are optimized through an outer loop. In this 
work, the use of the adjoint method for aeroelastic gradients enables simultaneous optimization of both sets of 
design variables, which is presumably more efficient. It should be recalled that, whatever control surface deflections 
are imposed by the optimizer, the aeroelastic analysis (Eqs. 1 and 2) will automatically locate the variables that trim 
the aircraft ( a and 6 for longitudinal maneuvers, a and /? for lateral). 

As seen in Figure 2, twenty control surfaces are attached to the trailing edge, and for the static aeroelastic 
analyses considered in this section, the attachments are rigid, and commanded flap rotations are exactly obtained 
(relative to the deforming wing). The same three load cases are used as in the previous section, but the third load 
case has both a symmetric (y 5 ) and an anti-symmetric ( y a ) piece, and their scheduling can be optimized separately. 
This results in 80 total DMCS design variables, in addition to the same structural design variables from above. Side 
limits of ±10° are imposed, and the scheduling design variables (as with structural design variables) are passed 
through a spatial filter in order to maintain a smooth deflection along the trailing edge. 

The optimal control surface scheduling is shown in Figure 8, where a positive flap deflection corresponds to a 
tip-down rotation of the surface. For maneuver load 1 (2.5g pull-up), additional lift is generated inboard, to 
maintain the overall lift-based trim of the aircraft, but load alleviation is used outboard, in order to reduce the total 
bending moment on the wing. This allows the optimizer to remove material from the structure (reducing mass, the 
ultimate objective function) without violating the strength and buckling constraints. Tip-up rotations of all of the 
control surfaces, from root to tip, would explicitly reduce the stresses and buckling loads even more, but the implicit 
effect would be an ultimate increase in a to maintain trim, and thus an increase in the aerodynamic loading. It is 
further noted that many of the 2.5g DMCS rotations reside at the side limits of ±10°; further optimization 
improvement could be expected if these limits were relaxed. 



*© — load case 1 : 2.5g 
-e — load case 2: - lg 
^ — load case 3: roll (sym.) 

-* — load case 3: roll (anti-sym.) 


Figure 8. Optimal DMCS rotations for each load case. 

Control surface scheduling for load case 2 (-lg push-over) has the expected opposite trends (but lower 
magnitude) as compared to the 2.5 g pull-up maneuver. For load case 3, the symmetric deflections are responsible 
(along with the angle of attack a) for maintaining steady level flight (N= 1), but these are largely set to zero by the 
optimizer. The anti-symmetric deflections are responsible for maintaining a steady roll rate (p ■ L/U = 0.065, as in 
all the cases in this paper), and so this is accomplished with a moderate tip-down rotation of all the control surfaces 
(of the right wing), with peak amplitude near the wing break. Inboard surface rotations are obviously less effective 
than outboard at generating a rolling moment, but they impose less torsional stress as well, and so the anti- 
symmetric scheduling shown in Figure 8 is the optimal compromise. 

The optimal structural design variables (upper and lower skin thickness and stiffener height) for this case are 
shown in Figure 9, along the span of the wing. This data is directly compared to the optimization case without 
DMCS design variables, which is repeated from Figure 4. The optimizer is able to significantly decrease the skin 
thickness, particularly in the upper skins (relative to the case without DMCS variables), though the reduction in 
stiffener height is milder. The mass of this wing is 4,242 kg, a substantial drop (31.3%) relative to the structure in 
Figure 4. 
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Figure 9. Optimal structural design variables (in the upper and lower skins) with and without DMCS design 

variables added to the optimization problem. 


VI. Open Loop Flutter Constraints 

The exercises in Sections IV and V have only included static aeroelastic considerations (stresses and buckling 
arising from trim loads), and will now be augmented with an open loop flutter constraint. This constraint stipulates 
that the wing not flutter within the flight envelope, including an additional margin based upon equivalent air speed. 
Flutter analysis is conducted using the same configuration seen in Figure 1, Figure 2, and Figure 3, though the 
aerodynamic paneling of the fuselage, tail, and engine are neglected. Furthermore, each control surface is assumed 
to be attached to the wing box with a flexible hinge spring. 

Aerodynamic analysis is computed using the doublet lattice method in the frequency domain. Open loop flutter 
constraints could also be computed in the frequency domain (using the p-k method), but are instead computed in the 
time domain (using the p-method). This is to provide a direct comparison with the results of the next section, where 
closed loop control requires state space modeling. The development of a state space aeroelastic flutter model is 
done with well-known methods [2], and will only be briefly described here. First, mode shapes (eigenvectors) of the 
flexible wing with rigidly-attached control surfaces are computed. Next, twenty control modes (one for each control 
surface) are computed as Ritz vectors, where each vector is zero except for a unit rotation of the corresponding 
control surface. The flexible mode shapes and the Ritz vectors are combined into a single modal matrix, and 
properly orthogonalized [25]. These modes are used to compute reduced mass, stiffness, and damping matrices, as 
well as generalized aerodynamic forces at tabulated reduced frequencies and Mach numbers. Frequency domain 
aerodynamics are converted to the time domain using the Roger approximation [26]. 

The final linear time invariant system is written as: 
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where b is half the mean aerodynamic chord, and M, K , and C are the reduced mass, stiffness, and damping 
matrices. A t are the generalized aerodynamic forces resulting from the rational function approximation, and y t are 
the corresponding lag roots (n total), x is the vector of modal amplitudes (wing deformation and control surface 
rotation), and x A is the vector of aerodynamic states. Eigenvalues of Eq. 5 are given by s = g + i - oo. For a given 
Mach number, these eigenvalues may be tracked across a range of matched point equivalent air speeds. Speeds at 
which a given eigenvalue crosses into the right-half plane ( g = 0) are flutter points. Any flutter point that occurs 
within the vicinity of the flight envelope by less than a 20% margin (as measured by equivalent air speed) is 
considered unacceptable. 
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The method of Ringertz [27] is used to formulate a flutter constraint, where the damping, g , of each mode is 
constrained to lie below a given curve at all equivalent air speeds U EAS of interest: 


0 0 < U EAS < U* 

S ■ ( U EAS - U*y U EAS > U* 


( 6 ) 


where 5 is a scaling parameter, and U* is the minimum allowable flutter equivalent air speed. Critical U EAS points 
(local minima) of the inequality in Eq. 6 are computed and lumped together into a single Kreisselmeier-Steinhauser 
constraint. As above, if this constraint is greater than one, then Eq. 6 is not satisfied, and the structure does not meet 
the required flutter margin. A separate KS constraint is utilized for each considered Mach number. Derivatives of 
the flutter constraint are easily computed with eigenvalue gradient techniques. A free-mode dynamic derivative 
approach is utilized here, where the derivative of the mode shapes with respect to design variables is neglected for 
the purposes of gradient computations. 

The optimization problem of the previous section is repeated here, adding the flutter constraint (at 0.85 Mach) 
to the static aeroelastic buckling and stress constraints. Both structural design variables and quasi-steady DMCS 
design variables are utilized as before. Several optimal results are presented in Figure 10 for various required flutter 
margins, and for two different DMCS hinge spring frequencies: 9 and 25 Hz. For reference, the first bending 
frequency of the optimal wing structures is on the order of 1 Hz. The dashed line in the figure indicates the 
minimum mass, which may be obtained with no flutter constraint: 4,242 kg, obtained in the previous section. As 
expected, enforcing a flutter constraint increases the optimal mass. For the 25 Hz hinge spring case (for which the 
control surfaces are nearly rigidly attached), the mass penalty is 34.4% under a 20% margin, which is typically 
required for civilian aircraft. For the more flexible hinge case, the penalty is nearly a factor of three. 



Figure 10. Pareto front between the flutter constraint boundary and the minimum structural mass. 

The aggressive flutter mechanism (and resulting mass penalty) associated with a flexible hinge spring is well 
known [28], and results from a situation where coupling between wing motion and control surface rotation 
exacerbates aerodynamic forces. More details are given in Figure 11, which tracks the various flutter mechanisms 
across a range of DMCS frequencies. This is done for the wing structure optimized without a flutter constraint. 
Three instabilities are seen for higher hinge frequencies, and the flutter speeds approach the values that would be 
obtained with a rigid hinge connection. These three flutter mechanisms become more severe (lower U EAS ) as the 
hinge spring frequency is decreased, and several additional mechanisms appear as well. The flutter speeds of these 
new mechanisms drop rapidly with decreased hinge frequency. At the lowest hinge frequency considered (7 Hz), 
the critical flutter speed is almost a third the value of the critical flutter speed computed with a rigid control surface 
assumption. 

Eigenvalue migration plots are shown in Figure 12 before and after the 20% flutter margin constraint is 
implemented, for a hinge spring frequency of 25 Hz. Similar data is given in Figure 13 for the lower value of 9 Hz. 
The top row of plots in both figures corresponds to the dashed line in Figure 10 (no flutter constraint), while the 
bottom row of plots corresponds to the data points with the largest mass in Figure 10. At the higher hinge frequency 
(Figure 12), the three flutter mechanisms noted in Figure 1 1 are seen here as well: two are hump modes, the other is 
a more severe instability. This is commonly seen for transport aircraft [29], where the primary hump mode is a 
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torsion mode dominated by engine pitch (mode 2), and the “hard” flutter point is a bend-torsion coupling. The 
dashed line in these figures is the constraint boundary of Eq. 6. This constraint is satisfied in the lower row of plots 
in Figure 12, and strongly active, as two modes brush up against the boundary in the vicinity of U*. 



Figure 11. Influence of control surface hinge frequency on the flutter mechanism. 




U EAS< m/S) U EA$ im/s) 

Figure 12. Eigenvalue migration without (top row) and with (bottom row) a 20% flutter margin constraint 

enforced: 25 Hz DMCS frequency. 

For the lower stiffness hinge spring (9 Hz, Figure 13), multiple hump mode and hard crossing flutter 
mechanisms are seen. The optimizer is also able to push each of these fluttering modes beyond the constraint 
boundary (lower row of plots), by adding a substantial amount of structure to the wing, as seen in Figure 10. Four 
separate modes are active in this case, but do not brush right up against the constraint boundary as seen in Figure 12. 
This is due to the conservatism of the aggregate KS constraint, which becomes more pronounced as more failure 
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mechanisms engage. Finally, the twenty 9 Hz modes (one for each DMCS) are also noted in the imaginary 
eigenvalue migration plots of Figure 13. 



Figure 13. Eigenvalue migration without (top row) and with (bottom row) a 20% flutter margin constraint 

enforced: 9 Hz DMCS frequency. 

Optimal structural design variables for the case with 25 Hz hinge springs and 20% imposed flutter margin are 
seen in Figure 14. Several key topological differences are noted here, as compared to the static aeroelastic optimum 
in Figure 4. Most predominantly, material is added to the wing tip leading edge, a clear attempt to alter the inertial 
loads via separation of coalescing modal frequencies. It is noted that landing loads and taxi bump loads are not 
included in the list of load cases here, but these types of inertial forces would penalize the tip mass [30] seen in 
Figure 14, and force the optimizer to be less reliant on this form of flutter mitigation. Secondly, the thickness of the 
spars is increased in the vicinity of the engine attachment points, which is presumably driven by the critical engine 
pitch fluttering hump mode. Finally, the optimizer is able to push most of the stiffener height design variables along 
the lower skin to the lower bound. The reason for this is unclear (large stiffeners are needed for the static aeroelastic 
case in Figure 4), though perhaps the added stiffness required to alleviate the flutter mechanism has implicitly 
decreased the need for large stiffeners as well. 

VII. Closed Loop Flutter Constraints 

It is seen in Section VI that an open loop flutter constraint can be satisfied using structural and quasi-steady 
DMCS design variables, though a weight penalty is imposed as well. For largely-rigid hinge spring connections (25 
Hz), the penalty is 34.4% under a 20% margin, and this penalty increases for more-flexible springs. This section 
explores the idea of relaxing the open loop flutter constraint, relying on closed loop control to meet the 20% margin, 
and assessing the weight savings from such a shift. At its most extreme, the open loop flutter constraint may be 
entirely removed, though this leads to a potentially dangerous situation where the vehicle relies on a controller to 
prevent a catastrophic aeroelastic instability within its flight envelope. More realistically (from a certification 
standpoint), the open loop margin may be relaxed to only 5% or 10%. Both cases will be considered here. Ideally, 
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the controller may be able to meet this flutter margin without needing to add extra stiffness (mass) to the wing, 
though the control cost (system power) may be excessive [2]. 



Figure 14. Optimal structural design variables for the optimization problem with static aeroelastic 

constraints and open loop flutter constraints. 

The aeroservoelastic linear time invariant system is given as: 


x ase — A ■ x ASE + B ■ u ASE (7) 

where x ASE is the state vector from Eq. 5, A is the matrix from Eq. 5, u ASE is a vector of input hinge moments (one 
per control surface), and B is a matrix which converts those hinge moments into generalized modal forces. The 
controller used here is a simple linear quadratic regulator (LQR), where feedback is assumed of the form: 

u ase = —K ■ x ASE (8) 


The feedback matrix is found which minimizes the performance index [31]: 


/ 



' Q ' X ASE + U ASE 


■ R ■ u ASE ) ■ dt 


( 9 ) 


where Q is the state weighting matrix, and R is the control weighting matrix. 

Two types of controller-centric design variables are considered here. The first are the members of the diagonal 
control weighting matrix R (one per control surface), which are passed through a spatial filter in order to maintain a 
smooth oscillatory shape along the trailing edge, much like the quasi-steady DMCS design variables. Upper and 
lower bounds are necessarily placed on the R design variables. Very large values (relative to Q, which is a fixed 
quantity for all cases in this paper) penalize the control action, and so instabilities will be arrested slowly. Very 
small values of R are indicative of cheap control, where instabilities can be damped out quickly. 

The second type of controller design variable is the flight condition (namely, the equivalent air speed U LQR ) at 
which the feedback matrix K is computed. Gain scheduling is the preferred long-term solution for active flutter 
control of transport wings [32], but for the conceptual studies provided here, a single control law is used for the 
entire flight envelope. The optimizer is able to choose the flight condition at which the controller is designed, 
however. A similar type of design variable is utilized in Refs. [6], [9], and [10]. Large values of the U LQR control 
design point (relative to the open loop flutter point) will stabilize the system at U LQR , as guaranteed by the LQR 
methodology [31], but the control cost / will be high. Small values of U LQR may not provide a closed loop flutter 
point measurably better than the open loop point, but the control cost will be low. 

The aeroservoelastic optimization problem solved in this section is to minimize mass, subject to the static 
aeroelastic stress and buckling constraints, an open loop flutter constraint, a closed loop flutter constraint, and a limit 
on the maximum allowable control cost /. The closed loop constraint is always set to a margin of 20%, but lower 
open loop margins are considered, as discussed above. Design variables include structural variables (shell thickness, 
stiffener thickness, stiffener height), quasi-steady DMCS rotation variables, the control weighting terms R , and the 
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controller design point U LQR . Gradients of the feedback matrix K and the performance index J are computed using 
methods in Ref. [33]. Once the K derivatives are found, derivatives of the closed loop flutter constraint are 
computed with the same eigenvalue sensitivity methods used for the open loop constraint. 

As with the open loop flutter constraint, a single Mach number of 0.85 is considered for the closed loop 
constraint. Furthermore, only the stiffer control surface hinge spring frequencies of 25 Hz are considered for this 
section. Prior to presented optimization results, trade studies demonstrating the dependence of the aeroservoelastic 
behavior upon R and U LQR are conducted. This is done for the static aeroelastic optimum found in Section V, with a 
mass of 4,242 kg. Figure 15 shows the closed loop flutter speed and the control cost for different values of R (the 
same fixed value for each control surface input). As indicated in the figure, the controller is designed at an unstable 
U L q R value 11.9% beyond the open loop flutter point. The LQR controller is guaranteed to stabilize the system at 
the design speed U LQR , even for very expensive controllers. For the highest values of R , the closed loop flutter 
speed is 35.2% higher than the open loop value. 



Figure 15. Dependence of closed loop flutter speed and LQR performance index upon weighting matrix R. 

For lower values of R in Figure 15, a higher premium is placed upon quickly arresting instabilities: substantially 
higher closed loop flutter points are seen in this area, nearly twice the open loop value. A caveat to this result 
(beyond the questionable idea of highly inexpensive control) is that some of the wind-off modes are unstable, as 
seen in the lower locus of flutter points in Figure 15. These modes re-stabilize (at a point not shown in the figure) as 
the equivalent air speed is increased to U LQR , as guaranteed by the LQR methodology. Similar behavior is noted in 
Ref. [9] and [34], and is presumably due to the stabilizing aerodynamic forces utilized by the controller at U LQR , but 
which are not present at very low dynamic pressures. A long-term solution to this problem is gain scheduling, 
which is not utilized here. Alternatively, the upper and lower bounds on R are (arbitrarily but necessarily) set at 50 
and 700, where the value of 50 is high enough that unstable wind-off modes do not appear during the optimization 
process. 

Similar data is provided in Figure 16 for changes in the LQR design condition U LQR , for fixed R values of 200. 
For the left plot, the system is always closed loop stable along the line with unit slope. For values of U LQR well 
beyond the open loop flutter point, the control cost / is very large and instabilities appear at lower equivalent air 
speeds, through a mechanism very similar to the previous figure. At very low values of U LQR , the closed loop flutter 
point approaches the open loop value, as the aeroelastic interactions that are used to derive the controller K are too 
weak to provide additional flutter mitigation beyond what is already available. 

Aeroservoelastic optimization results are now presented. Figure 17 shows the Pareto front between the optimal 
wing mass and the open loop flutter margin. The first curve shows a series of optimization results where only an 
open loop constraint is included (no control); this data is repeated from Figure 10. Next, the closed loop flutter 
constraint is set to 20%, and the open loop flutter constraint is progressively relaxed. This curve is coincident with 
the first, as the controller is able to meet the 20% closed loop constraint without adding any additional mass to the 
structure, and also without exceeding the performance index constraint limit, which is set to 10. For the final curve 
in Figure 17, the performance index limit is much harder to satisfy (set to one), and the LQR controller is unable to 
stabilize the system at a 20% margin without exceeding the allowable cost. The optimizer is then forced to add 
additional structure to the wing, with an 1 1.3% mass increase for the most critical case considered of 0.8 (-20% open 
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loop margin). Though not shown, the optimizer also re-adjusts the quasi-steady DMCS control surface rotations. 
These design variables have no explicit effect on either the wing mass objective or the flutter constraints, but can 
help maintain static aeroelastic constraints when material is re-allocated to satisfy dynamic constraints. 




Figure 16. Dependence of closed loop flutter speed and LQR performance index upon the LQR design point. 



* 0,L, flutter constraint only 

-e — C.L., performance index < 10 
— C.L , performance index < 1 


Figure 17. Pareto front between the open loop flutter constraint boundary and the minimum structural mass. 

Eigenvalue migration (real parts) for the -20% open loop margin case in Figure 17 is shown in Figure 18. Both 
open loop and closed loop eigenvalues are shown, and in both cases the corresponding flutter constraint (-20% 
margin for open loop, +20% for closed loop) boundary of Eq. 6 is strongly active. The closed loop case, in 
particular, has two modes closely interacting with the constraint boundary. Open and closed loop flutter modes 
(eigenvectors of Eq. 7) for this case are shown in Figure 19. The open loop flutter shape is a combination of wing 
bending, torsion, and engine pitch, and the control surfaces rigidly follow the trailing edge wing shape. The closed 
loop flutter shape shows less bending, and the input hinge moments (from the feedback controller) show noticeable 
control surface rotations relative to the wing, particularly in the region between the engine and the wing tip. 

The results of Figure 17 indicate a trade-off between structural mass and control cost, an idea discussed in Ref. 
[2] also. The complete trade-off curve between these two metrics is shown in Figure 20, where the open and closed 
loop flutter constraints are set to 0.8 and 1.2 (-20% and +20 margins), respectively. The corresponding control 
weighting design variables R along the trailing edge are shown in Figure 21. For higher allowable control cost 
values, the optimizer is able to meet the closed loop flutter constraint without adding additional mass to the wing, as 
noted in Figure 17 also. The R values for these cases (300) are unchanged from the initial value given to the 
optimizer, as these design variables are unneeded. Only when the allowable control cost constraint / is decreased 
below 2 is the optimizer finally forced to add mass to the wing, as this constraint is too restrictive for the FQR 
controller to satisfy on its own. 
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open loop dosed loop 




Figure 18. Open (left) and closed loop (right) eigenvalue damping under active flutter constraints. 



Figure 19. Open (top) and closed (bottom) loop flutter mode shapes. 

The expected monotonic increase in optimal wing mass with decreasing control cost requirements is seen in 
Figure 20 between / thresholds of 2 and 0.5. R values along the trailing edge are at the lower bound of 50, except 
near span stations that align with the engine attachment, where the upper bound of 700 is met. This is presumably 
related to the strong engine pitch noted in the closed loop flutter mode shape of Figure 19. Decreasing the cost 
constraint below 0.5 leads to strong discontinuities in the aeroservoelastic design space, whose source is ultimately 
unclear. There are no optimal designs between / thresholds of 0.5 and 0.1438, and again between 0.1438 and 0.084. 
The design at / = 0.1438 converges very slowly, and the optimization algorithm encounters discontinuities in the 
closed loop response and controller cost with small changes in the design inputs. This area of the design space is 
accompanied by a substantial change in the optimal R distribution along the trailing edge (Figure 21), where the 
inboard control surfaces are all at the upper bound of 700. A monotonic relationship between wing mass and / is 
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again seen below 0.1438, where every control weighting value R is pushed to the lower bound of 50. This is the 
expected value for a situation with very-stringent control cost limitations, given the near-proportional relationship 
between R and / seen in Figure 15. 



Figure 20. Trade-off between the allowable LQR control cost and the minimum structural mass. 



Figure 21. Optimal control weighting terms across the trade-off curve in Figure 20. 

VIII. Conclusions 

A series of aeroelastic optimization problems are solved in this paper, with increasing complexity in terms of 
constraints and design variables. All exercises are demonstrated on the uCRM wingbox, an un-deflected jig shape 
of the Common Research Model with conventional ribs, spars, skins, and smeared stiffeners. Panel method 
aerodynamics are coupled to shell-based finite element models of the built-up structure, and trimmed longitudinal 
and lateral maneuver loads are computed. Unsteady aeroelastic effects are computed in the time domain via a Roger 
transformation of the frequency- domain aerodynamics, and standard LQR controllers are utilized for closed loop 
mechanics. The following optimization problems are solved: 

1. Minimize wingbox mass subject to stress and panel buckling constraints spread across three static aeroelastic 
maneuver loads. Design variables include shell thickness, stiffener thickness, and stiffener height for design 
patches spread over the ribs, spars, and skins. 

2. Identical to (1), but structural design variables augmented with quasi-steady DMCS rotation variables, one 
per control surface (20) per static aeroelastic load case (3). Augmenting the design space with DMCS 
variables allows for an aeroelastically-feasible mass reduction of 31.3%. 

3. Identical to (2), but open loop flutter margin constraints included. This is done by ensuring that the real part 
of each state eigenvalue is below some constraint curve for all equivalent air speeds of interest, up to the 
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desired flutter margin. Satisfying a 20% flutter margin constraint increases the optimal mass by 34.4%, a 
penalty that can be much higher if the control surfaces are attached to the wingbox with flexible hinge 
springs. 

4. Identical to (3), but closed loop flutter margin and control cost (performance index) constraints included. 
Structural and DMCS design variables are augmented with control weighting design variables, as well as the 
equivalent air speed at which the LQR controller is designed. If higher control cost values are allowed, the 
LQR controller is able to meet the closed loop flutter constraint without the need to stiffen the structure. 
More stringent cost allowances do force the optimizer to add structural mass, however. 
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